Many-body approach to infinite non-periodic systems: 
application to the surface of semi-infinite jellium 
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A method to implement the many-body Green function formalism in the GW approximation for 
infinite non periodic systems is presented. It is suitable to treat systems of known "asymptotic" 
properties which enter as boundary conditions, while the effects of the lower symmetry are restricted 
to regions of finite volume. For example, it can be applied to surfaces or localized impurities. We 
illustrate the method with a study of the surface of semi-infinite jellium. We report the dielectric 
function, the effective potential and the electronic self-energy discussing the effects produced by the 
screening and by the charge density profile near the surface. 
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I. INTRODUCTION 

Density functional theory (DFT) is central to the 
understanding of several fundamental problems in con- 
densed matter physics, since it allows for obtaining the 
ground state properties of the system in a very accu- 
rate and computationally feasible way within an ab ini- 
tio framework. However, current approximations of the 
exchange correlation functional are based mainly on the 
Local Density Approximation (LDA), and hence on the 
properties of the homogeneous electron gas. Therefore 
they cannot account correctly for fundamental physical 
phenomena, like weak bonding, and consequently the van 
der Waals interaction. At surfaces, the LDA and derived 
Generalized Gradient Approximation (GGA) functionals 
are unsatisfactory and corrections for surface specific er- 
rors were recently proposed* 4 * Another example of the 
failure of the LDA/GGA is the potential calculated at 
large distances from a metal surface, which instead of 
reproducing the classical image one^ displays an expo- 
nential decay. This because the LDA/GGA neglects the 
long range electron correlations due to the inhomogcncity 
of the surface charge density. 

A way to improve on the limitations of the DFT solu- 
tion caused by approximated exchange-correlation func- 
tionals is to resort to methods of many-body perturba- 
tion theory. This also gives access to physical properties 
outside the realm of DFT, such as quasiparticle and col- 
lective excitations. 

The many-body problem was expressed by Hedin^ as 
a formally exact closed set of five equations that relate 
the single particle Green's function, the self-energy, the 
polarization, the effective two-particle potential and the 
vertex function. The GW approximation 4,5 (GWA) ne- 
glects vertex corrections and reduces Hedin's system to 
four closed equations. The results of DFT can be used to 
achieve an efficient scheme of solution of the GW equa- 



tions. In fact Dyson's equation for the exact Green func- 
tion requires the input of a reference one, which in the 
original Hedin's formulation is the Green function for the 
Hartree equation. A great improvement is to use instead 
the Green function of the Kohn-Sham (KS) equation, 
since it includes not only the Hartree potential but also 
exchange and correlation effects to some local degree of 
correctness. This minimizes the effort for self-consistency 
in the evaluation of the density. 

The GWA with the just mentioned input from DFT is 
a strategy that has been successfully applied to several 
systems. It combines the efficiency of the KS equation, 
with the use of the basic equations of many body theory, 
and it is able to supply reliable ground state properties 
as well as excitation energies, lifetimes, and generalized 
dielectric functions. For example the surface effective 
potential Vxc(z) with the correct image tail was achieved 
by Eguiluz et al.^ by solving the Sham-Schliiter equation 
in GWA, with a slab geometry. 

As already pointed out, DFT, sometimes coupled to 
the GWA especially for excited state properties^ have 
allowed for a realistic and accurate treatment of more 
and more complex systems of condensed matter physics. 
Among them, a single bulk impurity, solid surfaces, ad- 
sorbed molecules and clusters, interfaces as well as the 
recently investigated nano-contacts are indeed infinite or 
semi-infinite non-periodic electronic systems. However, 
in the theoretical description, such systems are usually 
confined in a fictitiously finite region/supercell with pe- 
riodic properties at the boundaries.- Consequently arti- 
ficial features may occur: for example, spurious inter- 
actions, non-physical oscillatory behaviors of the wave- 
functions, and a discrete spectrum in which it may be 
difficult to isolate localized electronic states and resolve 
resonant ones. For example, in order to properly account 
for the underlying bulk band structure, a DFT evaluation 
of the surface dielectric function for a semi-infinite crys- 
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tal, that avoids slab or supercell geometries, was recently 
proposed by Brodersen and Schattkei& 

In this paper we address the problem of the descrip- 
tion of an infinite non-periodic system within a GWA 
framework. We show here that the artificial reduction 
of an infinite volume to a periodic one can be avoided 
with minimal computational cost and more transparent 
results. Our method applies to systems (like surfaces or 
localized impurities) that asymptotically in space identify 
with others whose relevant correlators are well studied, 
like bulk crystal or vacuum. It relies on the assump- 
tion that, for any relevant correlator F(l,2), there is a 
finite space region Uf outside which F(l,2) cannot be 
distinguished from its asymptotic limit F°°(l, 2), within 
the required accuracy. This property is under control in 
the course of computation. It also provides a powerful 
tool for checking the correctness of the result, by veri- 
fying whether the correlator F matches with F°° at the 
boundaries of Uf or not. 

The virtue of this method is to replace the finite vol- 
ume of the cell or slab by a finite effective volume Up 
with boundary conditions that are intrinsic to the sys- 
tem. The property that "asymptotic" regions actually 
determine an effective region Uf of finite volume for the 
computation is connected with the principle of nearsight- 
edness recently introduced by KohnAiS It states that the 
local electronic structure near a point r, while requiring 
in principle the knowledge of the density (or the effec- 
tive potential) everywhere, is largely determined by the 
potential near r. 

In Sec.|n]we shall present the main points of the GWA 
to set the stage for further developments. Sec. lIIII outlines 
the method. Two main ingredients are the embedding 
approach for the zero-th order Green function, which 
guarantees that the properties of the infinite non-periodic 
system are taken correctly into account, and a Lemma for 
the inversion of infinite matrices. The application to the 
semi-infinite jellium is worked out in Sec. IIVI Because 
of its generality, the jellium surface is currently used as 
a bench mark system to evaluate many-body features. 11 
The extension of the GWA to the semi-infinite jellium can 
provide further data especially on how many-body prop- 
erties affect the spectral ones by calculating the dielectric 
function, the effective potential and the self-energy for a 
true continuum. Finally Sec.^is devoted to the conclu- 
sions. 



II. THE GW APPROXIMATION 

We wish to carry out a many-body treatment of an 
infinite non periodic system in the GWA. In order to 
present our work in a self-contained way, in this section 
we recall the basic properties and the equations of the 
GWA. To be specific, we write the Hamiltonian of a sys- 
tem of electrons with Coulomb interaction v(r,r'), in a 
static external potential V cx t(r) that couples to the den- 



sity n(r): 

k = E \& + E w ( f "^) + / drKxt(r)n(r). (1) 
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Atomic units (ao = 0.529 A, I hartree= 27.2 eV) are 
used throughout this paper. Since we are not interested 
in a spin-polarized phase, we consider a ground state with 
equal occupations for spin. The fermionic correlators are 
then proportional to the unit spin matrix. Two-point 
correlators, F(l,2) = F(riti,r 2 i 2 ) are time translation 
invariant and will be considered in frequency space: 



^(ri,r 2 , 



uj e 



-iu)(ti — 1 2 ) 



(2) 



The GWAASi is a self-consistent scheme that originates 
from truncating the exact closed set of five Hedin's equa- 
tions for the five basic quantities: Green function G, self- 
energy S, effective potential W, polarization P, and ver- 
tex function T. In the GWA the first two Hedin's equa- 
tions, namely the Dyson equations for the Green function 
G and for the the effective potential W, remain the same. 

The computational effort is reduced if, in the Dyson 
equation for G, one makes reference to the Green function 
Go that solves the KS equation: 
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Vks(r) 



Go(r,r» = <5(r-r'), 



(3) 



^ks (r ) = Vn (r ) + V sxt (r) + V xc (r) . (4) 

The Hartree potential Vn(r) = J dr'v(r, r')n(r') and 
the exchange-correlation potential Vxc( r ) contain the 
unknown density of the interacting system, which is 
to be found self-consistently by the relation n(r) = 
-2i f dujG (r,r,uj)e 1 ^. The factor e iur >, where tj -> 0+, 
results from time-ordering of operators and ensures ap- 
propriate analytic properties. 

It is simple to check that the exchange-correlation po- 
tential in Eq. Q modifies the Dyson equation for the 
Green function into: 

G(ri,r 2 ,cj) = G (ri,r 2 ,w) + J dr 3 J dr 4 G (ri, r 3 , w) 

x(£ X c(r 3 ,r 4 ,w) - Vxc(r 3 )<5(r 3 - r 4 ))G(r 4 , r 2 , uj). (5) 

Since the Hartree potential is accounted for exactly in the 
KS Eq. the self-energy diagrams are all of exchange- 
correlation type (with no tadpoles). The Dyson equation 
for the effective potential is: 

W(ri,r 2 ,uj) = u(ri,r 2 ) + 

y"dr 3 J dr 4 u(ri,r 3 )P(r 3 ,r 4 , w)VF(r 4 ,r 2 , w). (6) 

Next, there are the two Hedin's equations that basi- 
cally result from the unique skeleton diagrams^ of the 
exchange-correlation self-energy and of the polarization. 
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In the GWA they simplify, as one identifies the vertex 
with the bare one. 

£xc(ri,r 2 ,w) = 

f+OO J / 

■ / e ^iJ G (ri,r a) w W)W(r 2 ,ri,</), ( 7 ) 



P(Ti,T2,U)) = 

/+oo i / ^ 
— "G(n, r a> w + w')G(r 3 , r l5 a/).(8) 
-oo 27r 

The fifth Hedin's equation is neglected in GWA and 
corresponds to the skeleton expansion of the vertex, 
which is made of an infinite number of diagrams. The 
inclusion of vertex corrections is computationally expen- 
sive and difficult, with a degree of arbitrariness^ 

Finally we recall the following useful symmetry prop- 
erty of any correlator due to the invariance under time- 
reversal of the Hamiltonian: 

F{t 1 ,t 2 ,w)=F(t2,t 1 ,w) (9) 

For the polarization P and the effective potential W, 
Eq. implies further that they are even functions of 
the frequency u>. Such exact properties are preserved by 
the GWA. 



III. THE METHOD 

In this section we present a method to evaluate the 
interacting Green function of Eq. JSJ in the KS+GW 
scheme, for an infinite non-periodic system. In general 
in the GWA one has to perform the following operations: 

(i) Evaluation of the KS Green function Go in Eq. 

(ii) Evaluation of the ring diagram for Pq in Eq. JHJ. 

(iii) Evaluation of Wo by solving the Dyson Eq. J^J. 

(iv) Evaluation of S = iGqWq in Eq. J7J. (v) Evalu- 
ation of Gi by solving the Dyson Eq. JSJl . 

In principle, one should iterate the cycle (ii-v) for 
self-consistency, by inserting G\ in step (ii). Expe- 
rience with the homogeneous electron gas (HEG) has 
shown that spectral properties are better reproduced by a 
first-iteration calculation rather than by a self-consistent 

14 
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For periodic systems, the calculations can be restricted 
to a finite volume (e.g., the unit cell of a lattice crystal) 
with suitable boundary conditions. The present method 
implements the same procedure for infinite systems in a 
viable way. A weaker hypothesis is used: along the di- 
rections of broken symmetry any correlator F can be ap- 
proximated by the known F°° except for a finite length. 
Hence the calculation is carried on only in a finite volume, 
with boundary conditions determined by the "asymptotic 
properties of the system", i.e., by F°°. A direct space 
representation will be privileged in the directions where 
periodicity is absent. 



A. Preliminary steps 

We present two main ingredients which allow for the 
development of our method: the Green function embed- 
ding method to calculate Go in Eq. Q, and a Lemma 
for the inversion of infinite matrices, which permits to 
solve the Dyson equations for G\ and Wo in Eq. JSJ and 
Eq. JfJJ), respectively. 

1. The embedding method 

To start, the Green function Go of the KS equation © 
is needed. We make use of the Green function embed- 
ding methodi 15 i 16 i 17 Such a tool has been applied success- 
fully to the study of infinite systems without 3D period- 
icity, such as bulk impurities, surfaces and adsorbates>2> 
Its great advantage compared to the slab and the super- 
cell techniques is to provide a truly continuous density of 
states and the correct asymptotic behavior of all physical 
quantities close enough to the defect region. In the em- 
bedding method, space divides into a finite region V and 
one (or many) region V where the asymptotic regime is 
valid to the required accuracy. In this approach, the KS 
equation JS} in V U V is rewritten as an equation for the 
finite region V only. The effect of V appears as a sur- 
face term that adds to the KS potential. The modified 
KS equation, for r and r' in V reads as: 

[w-tf KS (r)]G (r,r»- 

/ d 2 r"[/ s (r,r'>)G (r",r>) = 

5(r-r'), (10) 

where Hks = — + Vks an d S is the boundary of V. 
The kernel Us(r, r', u>) is non-zero only for r, t'eS, where 
it is: 

C/ S (r,r» = i<5 s (r,r')n s (r) • V r + (Gg°)- 1 (r,r',c)(.ll) 

ns(r) is the unit vector normal to S in r, pointing out of 
V. G X3 (r, r', oj) is the Green function of the KS equation 
in the asymptotic region V', with Neumann boundary 
condition on S: 

[c-ff KS (r)]Gnr,r» - 6(t-t>), (12) 
n s (r)-V r Gg°(r,r',w) = 0. (13) 

(Gif 3 ) -1 is the surface functional inverse of G^, defined 
for points on S: 

[ dVG »(r 1 ,r', W )(GS o )- 1 (r',r 2 ,c) = < 5 s (r 1 ,r 2 ).(14) 

Eq. (|10f> can be expanded over a countable basis of V, 
thus reducing the evaluation of Go to a matrix inversion, 
which can be equally done for real or complex frequencies. 
We emphasize that the embedding method is formally 
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exact so that Go exhibits the truly continuous spectrum 
of the system. Being the solution of Eq. i|l(J|) , Go is known 
only for r and r' both in V. When the value of Go for 
one or both arguments outside V is required, it can be 
obtained with the "matching Green function" method^ 

2. Inversion of infinite matrices 

Consider the equations defining the inverse of a matrix 
A on two different volumes f2 and V, with V C O: 

/ dr 3 A(n,r 3 )yl f2 1 (r 3 ,r2) = <5(n - r 2 ), n,r a €0,(15) 
Jn 

/ dr 3 A(r u r 3 ) Ay 1 (r 3 ,r 2 ) = 6(n -r 2 ), n,r 2 € V(16) 
Jv 

In general Ay 1 is different from the restriction of Aq 1 
in V. However, the following Lemma gives a condition 
for the two matrices to coincide on a smaller subset U C 
V Cfl. 

Lemma: If (ri,r 2 ) = for all ri G fi — V and 
r 2 € U, then A^ 1 (ri,r 2 ) = A^ 1 (ri,r 2 ) for all ri,r 2 £ U. 
Proof: consider Eq. I|15l) for ri £ V and r2 G U, multiply 
it by Ay 1 (r4,r 1 ) and integrate in over V: 

/ dnAyXr^Ti) dr 3 A(r 1 ,r 3 )yl n 1 (r 3 ,r2) = A v 1 (r i ,r 2 ). 
Jv Jn 

The integral in r 3 over is split into an integral on V 
and on Q — V. The first integral yields A^ 1 (r 4 ,r 2 ), the 
latter vanishes because AQ 1 (ri,r 2 ) = for r\ G O — V 
and r 2 € U. • 

As a consequence of the Lemma, if we are interested in 
the values of A^ 1 in a subset U of the possibly infinite 
volume O it is sufficient to invert A on a suitable larger 
subset V, with V C O. Quite generally, the functions 
of interest have the property A^ 1 (ri,r2) — > as ri — 
r 2 | — > oo. Therefore, the hypothesis of the Lemma can 
be regarded as true to any degree of accuracy, for a large 
enough set V. 

B. Steps of the method 

1. The polarization 

In direct space and GWA, the polarization Pq is given 
by the ring diagram [Eq. (JHJ)] with Green functions Go- 

To compute it numerically, we note that the conver- 
gence of the integration in Eq. (JSJ is improved by us- 
ing AGo = Go — Gjf 3 , which decays faster than Go as 
| a; | — > oo. Hence we can take advantage of the known 
function Pq° , that corresponds to the asymptotic limit 
of Pq continued into the region of interest. We can write: 

P (w) = J dcu'[G^ + AG ](lu + uj') 

x[Gg° + AG ](^). (17) 



The integration of the term containing (Gj^Gg ) gives 
Pg 30 , the others have a faster decay for large \u>\. Owing 
to the non-analytic behaviour of the Green function close 
to the real axis, it is convenient to compute Eq. 1171 on 
the contour described in the Appendix. Since the spatial 
dependence is not involved in the computation of the 
polarization Po( r ij r 2,^), such dependence amounts to 
that of Go in the same region. 

2. The effective potential 

The effective (dressed) potential is the solution of the 
Dyson equation, Eq. ®, and it is formally given by the 
integral: 

Wo(ri,r2, w) = J dr 3 e _1 (ri,r 3 ,a;)w(r 3 ,r2). (18) 

where e _1 is the functional inverse of the dielectric func- 
tion e: 

e(ri,r 2 ,w) = <5(ri - r 2 ) 

-J dr 3 u(ri,r 3 )Po(r 3 ,r 2 ,w). (19) 

The decay properties of the polarization Po as ri — r 2 \ — > 
oo imply that Eq. (|19fl can be evaluated numerically by 
introducing cutoffs for the integration variable r 3 . 
The inverse dielectric function is formally defined by: 

J dr 3 e(ri,r 3 ,w)£ _1 (r 3 ,r2,w) = S(r 1 - r 2 ). (20) 

The inversion of the matrix e on an unbounded region is 
not a feasible numerical calculation. However, we only 
need evaluate e _1 on a finite region U £ outside which the 
asymptotic regime holds. Since in ordinary systems the 
effective interaction between electrons far apart decays 
to zero as the distance increases, Eq. (|18|) implies that 
also e _1 (ri,r 2 ) goes to zero as |ri — r 2 | —* oo. So we can 
use the Lemma proven in Sec. IIII A 21 and restrict the 
integration in Eq. I|20l) to a finite region V £ , U £ C V e , still 
obtaining correct values of e _1 (ri,r2) for ri,r 2 G U e . 
An expansion over a discrete basis set is now possible, 
leading the problem to an ordinary matrix inversion. The 
size of the region V e is simply a numerical parameter, and 
convergence in the resulting must be checked. With 
this purpose, a localized basis set [like 4>i(r) = c5(r — r^)] 
is more convenient. In this case the matrix elements , 
being directly proportional to e -1 ^,!^), do not change 
except for the normalization factor c, when the size of V s 
or the number of basis function changes. The effective 
potential can now be evaluated from Eq. (|18fl . We only 
remark that the frequency argument to is fixed, so that 
Wo is computed for the same values of w as those chosen 
for the polarization - u> purely imaginary as shown in the 
Appendix. 

Finally we note that the assertion that e _1 (ri, r 2 ) — > 
as |ri — r 2 | — > oo before e -1 is actually evaluated does 
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not pose conceptual difficulties: the decay to zero has in 
fact to be checked for the known function (e°°) _1 . 



3. The self-energy 

It is customary to split the self-energy into the sum of 
the exchange Sx = iGv and the correlation term Y,q = 
iG(W — v). The evaluation of the exchange term poses 
no problem: 

S x (r 1) r 2 ) = w(r 1 ,r 2 ) J °° ^- e ^G Q {v u r 2 , w). (21) 

The frequency integration sums over the occupied KS 
states, and can be performed analytically if the spec- 
trum is discrete. In our case, the spectrum is generally 
continuous [Go is the solution of Ea. i(TU|l ]. and the sum 
is replaced by an integral to be performed numerically. 
A change of contour proves useful, together with the in- 
formation that we shall be dealing with systems in which 
there are no KS states below the bottom of the band Eq 
(i]' is a positive infinitesimal): 

/+oo 
dwe M "GoH 
-oo 

= f 6u[G (u-ir/) -G (w + *»/)]• (22) 

The resulting integral is computed straightforwardly, /i 
is the KS chemical potential. 
The correlation term is: 

/+oo 
dct/e iaA? G (w + to') [Wo(w') - v] . (23) 
-oo 

Since the main contribution to the above integral is 
Eg 3 , one can use the same splitting of the self-energy 
in Eq. Il2.'ill into a leading asymptotic term plus a correc- 
tion due to the inhomogeneity as in Sec. 1111 B II for the 
polarization. 

Mathematical details on the calculation of the inte- 
grals determining the the self-energy are discussed in the 
Appendix. 



4- The Green function 

The Dyson equation for G [Eq. JSJ] is formally the same 
as the Dyson equation for W [Eq. once G, Exc — Vxc 
and Go are identified with W , P, and v respectively. 
Therefore we define the function exc( r i, r 2, w), which is 
analogous to the dielectric function: 

exc(ri,r 2 ,w) = 6(rx - r 2 ) - J dr 3 Go(rx, r 3) w) 
x[E X c(r 3 ,r 2 ,w) - V5cc(r 2 )<5(r 3 - r 2 )] . (24) 



Even with the these strong analogies, exc presents a 
striking difference with e: the decay of exc( r i, r 2) as 
|ri , r 2 1 goes to infinity is linked to the one of Go , which in 
turn varies according to the dimensionality of the system. 
As a consequence, exc( r ij 1*2) may not go to zero as |ri — 
r 2 1 goes to infinity, thus not satisfying the hypothesis of 
the Lemma in Sec. IIII A 21 To make Go decay as |r — 
r'| — > 00 (which implies the same property for exc)> h 
is convenient to solve Eq. (|2~4*}l at a complex frequency 
lu + iA. The choice of the real quantity A depends on 
a compromise: if it is too small, the decay of Go is very 
slow and a large region of inversion of exc is needed; if 
it is too large, the structures on the real frequency axis 
we are interested in are broadened - in fact A plays the 
role of resolution. The final result of the calculation, 
the interacting Green function, is thus evaluated on a 
translated frequency axis uj+iA. Analytical continuation 
improves the resolution: first, the values of G are fitted 
with a rational function, then the expression is continued 
to the real frequency axis. 



IV. SEMI-INFINITE JELLIUM 

In this section we illustrate the application of the 
method to semi- infinite jcllium. 19 



A. Introduction to the system 

Semi-infinite jellium is a neutral system of dynamical 
electrons in a background of uniform positive charge den- 
sity in the half-space z < 0. In the half-space z > there 
is no positive charge. We denote the uniform density by 
n, n = l/(47!Tg/3). In the computations the value of r 8 
will be fixed to give an electron density equal to that of 
aluminium (r s — 2.07 ag). 

Besides time translations and time reversal, the system 
is also invariant under translations parallel to the surface. 
Hence the wave- vector parallel to the surface k|| is a good 
quantum number. It is important to observe that since 
the solid is semi-infinite, the wave-vector k z may take 
any real value. So we shall deal with truly continuous 
densities of states, also for a fixed ky. A slab calcula- 
tion would instead determine a discrete spectrum (i.e., a 
nonphysical quantized set of k z wave- vectors) . Because 
of the just mentioned symmetry, a two-point correlator 
F(r|| , z, ; rj, , z'; u>) only depends on the difference r\\ — rj. , 
where rii = (x,y). We can express F in terms of its 
Fourier transform (FT) with respect to a wave vector 
parallel to the surface: 

F(r {{ t z; rf,, z'; u ) =j S e<k "' ( '" ~^ F ^ k|| , w).(25) 

Since the system is also invariant under rotations around 
the z-axis, F depends only on the modulus of ky . 
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The study of semi-infinite jellium is basically one- 
dimensional and satisfies the requirements of applicabil- 
ity of our method. In fact, the perturbation induced 
by the jellium edge is localized near the surface. At a 
distance of few r s inside the solid, the properties of the 
system approach those of the infinite, homogeneous elec- 
tron gas (HEG). Many-body results for the HEG in the 
GW approximation are well known in the literature^ If 
we indicate by F™ G (fc, uS) the correlator F evaluated for 
a HEG of density n, the bulk and vacuum limits of F are 

F^ y (z, z', , W )= J ^MW) F HEG ( ^ k | + ^ ul2Q) 

Therefore we need evaluate the correlator F only on a 
limited interval [2b,2v] of the z-axis, which will be the 
volume Uf of the previous Sec. IIII A 21 

Care must be taken when F is evaluated at z and z' 
in different regions of space, e.g., z in bulk and z' in 
vacuum. However, we are working with functions with 
the important property F(z, z') — > when \z — z'\ — > 00. 

I 



where v is the FT of Coulomb's potential: 

2-7T 

v(\z 1 -z 2 \,k ]] ) = -^e- k ^-^. (28) 

fc|l 



C. Steps of the calculation and results 

We omit details regarding the starting point: the KS 
Green function Gq(z, z' , ku , u>) can be obtained straight- 
forwardly, with the embedding method described in 
Sec. IIII A II We calculated it in the LDA. 



This guarantees that when z is in bulk and z' in vacuum 
the functions Fb(z, z') and Fy(z, z') are both zero to the 
desired accuracy, provided that |zb — Z\r\ is large enough. 

The dependence of F on its four arguments (z, z', kn , cj) 
is in principle continuous. Numerically, we stored the in- 
formation about F on a four dimensional discrete mesh. 
Intermediate values, when necessary, are obtained with 
interpolation algorithms. For z and z' meshes, natural 
limit values are given by zb and zy. Cutoffs for fey and uj 
can be fixed since F — ► when fey — > 00 or \u>\ — » 00. Dif- 
ferent meshes have to be chosen for different functions. 
Numerical convergence must be checked for any param- 
eter defining the mesh. 



B. G0W0 equations 

In the chosen representation, the Go Wo equations 15181) 
take the form: 



(27a) 



(27b) 



(27c) 



(27d) 



1. The polarization 



The parallel-wavector convolution in Eq. I|27a|) does 
not pose numerical difficulties. Regarding the frequency 
convolution, the factor e lu} v is necessary for the conver- 
gence of the integral. In fact, the Green function Go(uj) 
approaches zero like |w| _1 / 2 as \uj\ — > 00 (this can be 
easily verified for the HEG Green function). We follow 
the treatment of Sec. IIII B II and separate the integrand 
in two terms: one (G™ G G™ G ), that gives P™ G , and 
(G G - G™ G G™ G ) that can be evaluated numerically, 
because it goes to zero as |w|~ 3 / 2 when \uj\ — > 00. 



r +oc dcu' r d 2 k(, . , 

Pb(zi,z 2 ,k h u) = -2i J —J j^e™«Go{z 1 ,Z2AH+k\\\,u + oj')G (z2,zi,k' p u ) '), 



W (z 1 ,z 2 ,k\\,uj) = v(\zi - z 2 \,k\\) + dz 3 / d^d^i - z 3 \,k\\)Po(z3,Z4, k\uui)Wo(z4,,Z2,kn,u)), 



+00 dLJ > r d 2 k', 



Sxo(«i,«a,*||,w)=t J —J j^e^' 1 G (z 1 ,z 2 ,\k ll +-k\ l \,u + uj')Wo(z 2 ,z 1 ,k\ l ,Lu'), 



G 1 {z 1 ,z 2 ,k\\,w) = G (zi,Z2,k\\,u}) + J dz 3 j dziG {zi,z 3 , fey, w) 

x(Y>Kc{z?,,z 4l k\\,u)) - Vxc(z3)S(z 3 - z 4 ))Gi(z 4 , z 2 , fey , uj), 
I 
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2. The effective potential 

As shown in Sec. IIII B 21 one has to evaluate the di- 
electric function e, defined by 

e(zi,Z2,k\\,w) = 5(zi - z 2 ) (29) 

-J dz 3 v(\z! - Z3\,k\\)P Q (z3,Z2,k\\,u), 

and its inverse e _1 in order to solve Dyson equation for 
W, Eq. (l27bt : 

W / (zi,z 2 ,fc||,w) = 

y dz 3 £~ 1 (zi,2 3 ,fc||,Lj)u(|z 3 - Z 2 |,/C||). (30) 

The z integrations may be restricted to finite intervals 
since Pq(zi,z 2 ) and e~ 1 (zi,z 2 ) go to as \z\ — z%\ — > oo. 
We also point out that the calculation time required by 
the integrals in Eqs. (|29|l and (|30[l is negligible, so the 
cutoffs of z can be cautiously overestimated. 

To obtain e -1 , defined by the inversion over the whole 
z-axis, 

/+oo 
dz 3 e^ 1 (z 1 , z 3 )e(z 3 , z 2 ) = S(z 1 - z 2 ), (31) 
-oo 

on the finite interval U e — [zb, zy], we exploit the Lemma 
in Sec. IIII A~2l and restrict the integration in Eq. j|31|l to 
the finite interval V e — [Lb,Lv], U e cV e . 

Correct values of e~ 1 {zi 1 z 2 ) for z i, z 2 G U £ are ob- 
tained if Lb and Ly are conveniently large. It is suffi- 
cient that e" 1 (z, zb) — if z < Lb and e^ 1 (z, Zy) ~ if 
z > Ly to the required degree of accuracy. This concept 
is graphically presented in Fig. ^ where we show e -1 for 
different values of Lb and Ly. The negative peak for 
Z\,Z2 — Lb (on the left of each plot) represents a spuri- 
ous feature introduced when the region of integration is 
restricted to a finite interval. Notice that this behavior 
is located at the boundaries of V e regardless to its size. 
Since e~ 1 (zi, z 2 ) is different from zero only for values of 
z 2 close to z\ , this nonphysical feature is confined within 
few atomic units from the boundaries. Hence the inter- 
val V e has to be only slightly larger than U s . For the 
values described in Fig. ^ if z b = — 15 do, the choice 
Lb = —20 ao is already an accurate one. A similar dis- 
cussion has to be done with respect to Ly, but in this 
case the spurious peak is much smaller. 

In Fig. |21 we show the contour levels of the difference 
between the effective and bare interaction Wq — v in the 
Z\,Z2 plane. The HEG levels of the same quantity are 
also reported. The agreement is excellent when z% and 
z 2 approach bulk. As we move into the vacuum, Wq — v 
correctly goes to zero. 

Next we consider the effective potential Wo in the more 
intuitive representation, (ri,r 2 ,uj), obtained by anti-FT 
with respect to k|| . For simplicity, we limit our discussion 
to the static case (u> = 0) and consider collinear points 



fei. z 2 ) ~ 5 (z r z 2 ) 




FIG. 1: Values of e~ 1 (zi, Z2) — b~{z\—Z2) for different inverting 
regions V e = [Lb,Lv]- Here, r s = 2.07 ao, w — and k\\ = 
0.2kF, kF being the Fermi wave-vector. 




-10 -5 5 



Distance from the surface Z\ {%) 

FIG. 2: Contour levels of WoO^i: Z2, k\\, u>) — v{\z\ — zi |, fcy) 
in the 21,2:2 plane. Thick curves: semi-infinite jellium. Thin 
curves: HEG of equal density. Here r„ = 2.07 ao, w = and 
fc|| = 0.2fe F . 



on the normal to the surface (r^i = r 2||)- Fig. |21shows 
the effective potential as a function of z from bulk to 
vacuum: Wq is similar to a Yukawa screened potential for 
z\ and z 2 in bulk, and it coincides with the bare Coulomb 
interaction for z\ and z 2 in vacuum. Some intermediate 
values are shown: for Z\ fixed near the surface, Wq is 
no longer a symmetric function of z 2 with respect to z\, 
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o 

tj 




-10 -5 5 10 

Distance from the surface z 2 (.%) 

FIG. 3: Effective and bare Coulomb interactions near the 
jellium surface for points aligned on the normal. r s = 2.07 ao- 



being the screening inhomogeneous. 



3. The self-energy 

The evaluation of the self-energy follows closely that of 
the polarization. However, differently from the previous 
case, it is not necessary to split the self-energy Ec into 
a leading term and a correction to it, as discussed in 
Sec. IIII B 31 because the difference W — v decays fast 
enough to ensure convergence. 

In Fig. 0] we report the contour levels of the self-energy 
evaluated at ku = and u> = fi . Fig.^Ja) shows the result 
of our calculation: a particular feature in the near-surface 
region is the "arabian" shape of the contours levels. This 
arises because when z% and z 2 lie outside jellium, the self- 
energy, as \z\ — Zi\ increases, decreases in a slower way 
than in bulk owing to lower screening. 

In order to reduce the numerical effort required by a 
GW calculation, attempts have been made to mimic the 
self-energy with efficient models. In Fig. 0{b) we report 
the self-energy of a model based on fitting the spatial 
dependence of £ from that of the HEG having the av- 
erage density of the system, successfully tested for bulk 
siliconiSi The surface peculiar contours in Fig. 0Ja) can- 
not be reproduced by this model and the decay of £ in 
vacuum is much slower. The comparison between these 
two results points out that a model based on an average 
density cannot reproduce features which are characteris- 
tic of a surface. 

In Fig. 0Jc) we show the self-energy evaluated from 
a Go Wo calculation where, instead of the self-consistent 
DFT-LDA Green function Go, one uses the Green func- 
tion Gg tep computed for the approximated effective po- 
tential: vfff P (z) = vxc{— 00 ) if z < and v s c t g P (z) = 
otherwise. It is well know that in this case there are 
stronger Friedel's oscillations, and the electronic states 



-5 




0.025 
0.050 
0.100 
0.200 



-5 5 

Distance from the surface Zj (a ) 




FIG. 4: Self-energy Exc(«i, «2, w) for r s = 2.07 a , k\\ = 
and u) = fi. (a) Result of the GW computation, (b) Paula 
Sanchez-Friera's model, (c) Step-potential approximation. 
Contour levels as in (a). 



(and thus the Green function) decay faster away from the 
surface. This is reflected in the self-energy. Finally we re- 
mark that the difference between G Wo and G tep W step 
is due primarily to the difference between the Green func- 
tions, whereas the effective interactions are very similar. 



4- The interacting Green function 

Following Eq. I|24|l . the interacting Green function Gi 
for the semi-infinite jellium can be evaluated from 

Gi{z u z 2 ,k\\,uj) = 

dz 3 e x ^(zi, z 3 , k\\ , uj)Gq(z 3 , z 2 , k\\ , w), (32) 

where the kernel e^c is obtained by inverting 



exc(zi,Z2, k\\,w) = 5(zi - z 2 ) - / dz 3 G (z ll z 3 ,k\\,uj) 
x[Exc(z3,z 2 ,k\\,oj) - Vkc (22)5(23 - 22)]- (33) 
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When the frequency is infinitesimally close to the real 
axis and its real part corresponds to extended states, 
Gq(zi, Z3, feii, lu), as a function of Z3, is a plane wave 
with undamped oscillations propagating to infinity. As a 
consequence e^ c (zi, z%) does not go to zero as \z\ — Z2I 
goes to infinity, thus not satisfying the hypothesis of the 
Lemma. One must then recur to the procedure outlined 
in Sec. 11111341 i.e., solve Eq. Ij27d(l at a complex frequency 
lu + iA. We experienced that an interval V e about 100 a 
wide was needed for a value of A of about 0.05 hartree, 
in order to describe the surface region correctly. 

The imaginary part of the interacting Green function 
G yields the many-body spectral weight function 

A(z, kii , lu) = — — SsG(z, z, kii, w)sgn(aj — ^). (34) 

11 ^ 

This quantity gives a measure of the quasi-particle am- 
plitude and is directly related to a variety of experi- 
ments such as photoemission spectroscopies^ and scan- 
ning tunneling microscopy^ The integral in k|| gives the 
local density of states (LDOS) 



a(z, lu) = 




The evaluation of the LDOS of semi-infinite jellium 
in this framework demonstrates the presence at the sur- 
face of a broad image-potential induced (IPI) resonance, 
which emerges sharply when results are compared to 
DFT-LDA ones. We stress that an IPI resonance width 
can only be obtained by a many-body approach like ours 
which takes into account the semi-infinite character of the 
solid. We refer to Ref. |2j| for the results and a detailed 
discussion on this topic. 

V. CONCLUSIONS 

We have presented a method to investigate infinite 
non-periodic systems in the framework of the GWA. Cal- 
culations can be performed in finite regions, without in- 
troducing nonphysical boundary conditions, such as con- 
fining barriers (the slab approach) or a 3D fictitious pe- 
riodicity (the supercell one). In such systems (e.g., a 
solid with a surface) densities of states are continuous, 
and while really discrete states may exist inside gaps, 
other ones become resonances when they do overlap in 
energy with a continuum band. The proposed method 
is particularly suitable for the description of these sys- 
tems. In fact on the one side the embedding approach, 
which allows for calculating a truly continuous density of 
states, includes automatically the hybridization between 
bulk and surface states. On the other many-body effects, 
whose treatment is needed for excited states or image po- 
tential ones, are accounted for at the GWA level. 

On the contrary a DFT slab calculation of such sys- 
tems (e.g., in the LDA or GGA) is only able to work out 
a spectral weight constituted by delta functions, one for 



each discrete eigenstate, while the real structure of the 
spectrum may be in general more complicate as just out- 
lined. The GWA correction cannot amend by itself this 
result, but only determine a broadening of quasi-particle 
states (plus eventually minor additional structures) due 
to many-body correlations. This broadening, which can 
be evaluated in first approximation by taking the aver- 
age value of the self-energy over the DFT state, may be 
much smaller than that due to hybridization effects, as 
it is the case for IPI resonances. 

In this paper we have also extensively investigated 
semi-infinite jellium by our approach. We have illus- 
trated the bulk-to-vacuum transition of the many-body 
electron gas properties. By comparing the LDA and 
GWA density of states, this method has been able to 
identify an image potential surface resonance of large 
width. 2 - Extension of this approach to semi- infinite real- 
istic surfaces could bring a wealth of accurate data on the 
spectral properties of surfaces and adsorbates, especially 
regarding the excited states. 
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APPENDIX: ANALYTIC CONTINUATION OF 
FREQUENCY INTEGRALS. 

The presence of non-analyticities close to the contour 
of the frequency integration renders it difficult to inte- 
grate expressions containing the Green function G and 
the effective potential W numerically, as for the polar- 
ization or the self-energy. 

Consider the integral defining the polarization P in 
Eq. (JSJ first. The Green function has poles (or cuts) 
just below the real lu axis for u> > \x and above for 
lu < fi. Therefore, if z is a pole or a point in the cut, 
sgn(/i — Kz) = sgn(3z). Note that the factor e lu v im- 
plies that only the residues related to occupied states 
(lu < /i) axe summed up. To avoid the just mentioned 
numerical difficulty one can define the analytic contin- 
uation of P to complex frequencies as the sum over the 
same residues, now evaluated at a complex frequency »S4 
It is easy to show that for purely imaginary frequencies 
this corresponds to rotate the integration contour to the 
complex frequency axis fj, + iu' (11,11' real). In the GWA 
the continued polarization [Eq. (JSJ)] is: 

P (iu) = -2i do/G (w' + iu)G (u'). (A.l) 

J fi—ioo 

On the same footing, also the self-energy [Eq. Q] can 
be continued to complex frequencies. If lu = fj, + iu, the 
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following relation holds: 

/ioo 
duj'G {uj' + v + iu)W Q {uj')(K.2) 
-ioo 

where the analytic continuation of Wq{uj' — iu') is eval- 
uated by inserting Po(iu) into Dyson's equation. Note 
that the Lchmann representation of Pq implies that a 
pole z of Wq has sgn(3z) = — sgn(5ftz). 

The self-energy resulting from Eq. (|A.2J) will be known 
on the complex line to = fi+iu. This is useful for the eval- 
uation of integral properties (e.g., the total energy), but 
for spectral properties the Green function (and hence the 
self-energy) has to be known at real frequencies. To this 
end, one can fit £xc on the complex axis with a simple 
analytic expression, to be continued to the real axis^ 
The multi-pole expansion is perhaps the more common 
one: 



N 



SxcM = fl o + ^2 ■ 



(A.3) 



A small number of poles (N = 2^4) normally provides 
a good fit. 

Finally we recall a useful result to rotate the integra- 
tion path in frequency space. Consider the two integrals, 
where lu, a and b are real: 



Fi(w) 



dw'- 



1 



(a/ - zi)(u + lu' - z 2 ) 
sgn(3?zi) - sgn(3,Z2) 



Z2 + 21 



(A.4) 



a+ioo 



dw'- 



1 



(u)' — z\){b + iuj + uj' — z 2 ) 
sgn(a — $tzi) — sgn(a + b — %lz 2 



b + iui — z 2 + z\ 



(A.5) 



The two numerators are equal if sgn^zi) = sgn(a — 
$tzi) and sgn(3z2) = sgn(a + b — §lz 2 ). In this case: 
F 2 (tu) = F\{b + iuj), i.e., F 2 is the analytic continuation 
of Fi to complex frequencies b + iu>. Notice that, to be 
analytic, the continuation has to be performed after the 
integration. 



If both z\ and z 2 are poles of the time-ordered Green 
function [as for the polarization, Eq. ((SJ], it follows from 
the Lehmann representation that the condition is met for 
a = /i and 6 = 0. If z\ is a pole of the effective interaction 
and z 2 is a pole of the time-ordered Green function [as for 
the self-energy, Eq. JJJ], the condition is met for a = 0, 
b = fi. 
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